Introduction to Python Day 5

Verjinia Metodieva and Daniel Parthier

2025-03-18

SciPy

A python library for scientific questions
NumPy in the background
Many useful functions for signals (filtering, finding peaks, FFT, deconvolution)
Some statistical tests (t-test, wilcoxon, correlation, etc.)

Signal processing

Module signal
Important part of analysis
Getting information of interest

Detect peaks

Function to detect peaks in data: find_peaks
You can specify criteria
- Width (minimum, maximum)
- Distance to next peak
- Threshold
- etc.

Detect peaks

import json
import numpy as np
import matplotlib.pyplot as plt

sampling_rate = 2e4
with open('data/charact_data.json') as f:
    data = json.load(f)
# print the keys
data = {key: np.array(data[key]) for key in data.keys()}
time_s = np.linspace(0, data['D1'].shape[1]/sampling_rate, data['D1'].shape[1])

Detect peaks

plt.plot(time_s, data['D1'][10])
plt.show()

Detect peaks

from scipy import signal
peak_location, peak_properties = signal.find_peaks(data['D1'][10], height=0, width=[10, 50], distance=10)
print(peak_location)

[ 3310  3927  5223  6880  8638 10530 12395 14431 16460 18490 20600]

Detect peaks

plt.plot(time_s, data['D1'][10])
plt.plot(time_s[peak_location], data['D1'][10][peak_location], '.')
plt.show()

Detect peaks

Peak properties:
- Peak height, widths, prominences, thresholds
- Depends on function input
- More details in documentation

Practice

Find the peaks of all traces
Append the location and height to separate lists
Plot the number of spikes against the current
- [-300, -200, -150, -100, -50, 0, 50, 100, 150, 200, 250, 300]
Think of a good way to make use of functions

Fit data

Modul optimize
Different functions to fit lines/functions
curve_fit allows to specify a function which is fit to data

Specify a function

def func(x, a, b, c):
    return a * np.exp(-b * x) + c

x = np.arange(0,10,1)
y = func(x, 1,0.3,0.1)
plt.plot(x, y, 'o')
plt.show()

Fit function

from scipy import optimize
popt, pcov = optimize.curve_fit(func, x, y, p0=(1,1,1))
x_fit = np.arange(0,10,0.1)
plt.plot(x, y, 'o')
plt.plot(x_fit, func(x_fit, *popt))
plt.show()
print(popt)

[1.  0.3 0.1]

Convolve and Deconvolve

Signals can be distorted
Multiplying something to signal with convolve
Remove known distortion from signal deconvolve
Where do you see use cases?

Signal analysis (power)

time = np.linspace(start=0, stop=0.5, num=2000)
mu, sigma = 0.25, 0.01
sinewave = np.sin(time * 250 * np.pi)
gaussian = (1 / (np.sqrt(2 * np.pi * np.square(sigma))) *
            np.exp(-(np.square(time - mu) /np.square(2 * sigma))))
ripple = gaussian * sinewave

Signal analysis (power)

np.random.seed(0)
trace = ripple + gaussian*6 + np.random.normal(0, 3, size=time.shape)
plt.plot(time, trace)
plt.show()

Signal analysis (power)

Module signal
Power spectrum to analyse frequencies
e.g. welch

Power spectrum

f, Pxx_den = signal.welch(trace, fs=1/time[1], nperseg=1024, average='median')
plt.semilogy(f, Pxx_den)
plt.vlines(125, Pxx_den.min(), Pxx_den.max(),color='r', linestyle='--')
plt.xlabel('Frequency [Hz]')
plt.ylabel('PSD [V**2/Hz]')
plt.title('Power Spectral Density')
plt.show()

Returns frequencies f and power densities Pxx_den

Power spectrum

Text(0.5, 0, 'Frequency [Hz]')

Text(0, 0.5, 'PSD [V**2/Hz]')

Text(0.5, 1.0, 'Power Spectral Density')

Filtering data

Filtering data is used to remove different components from signal
lowpass, highpass, bandpass
different types of filters (Butterworth, Chebyshev, FIR, etc.)

Filter construction

Need different parts
- Nyquist frequency
- Filter frequencies (0 to 1 - frequency/Nyquist)
- Type of filter
- Order (‘strength’)
Outputs are filter coefficients

from scipy import signal
sampling_rate = 1/time[1]
nyquist = sampling_rate / 2
b, a = signal.butter(N=5, Wn=[100/nyquist, 250/nyquist], btype='band')

Filter construction

w, h = signal.freqz(b, a, fs=sampling_rate, worN=2000)
plt.plot(w, abs(h))
plt.title('Filter frequency response')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Gain')
plt.grid()
plt.axvline(100, color='r', linestyle='--')
plt.axvline(250, color='r', linestyle='--')
plt.show()

Filter construction

Text(0.5, 1.0, 'Filter frequency response')

Text(0.5, 0, 'Frequency [Hz]')

Text(0, 0.5, 'Gain')

Apply filter

Different ways of filtering (IIR/FIR)
one way or two ways (filtfilt)

ripple_filtered_trace = signal.filtfilt(b, a, trace)
plt.plot(time, ripple_filtered_trace)
plt.plot(time, ripple, '--')
plt.show()

Filter for lowpass

b, a = signal.butter(5, 80/nyquist, btype='low')
lowpass_filtered_trace = signal.filtfilt(b, a, trace)
plt.plot(time, trace)
plt.plot(time, lowpass_filtered_trace, '--')
plt.show()